## Preamble #### 
library(here)
library(plm)
library(lmtest)
library(sandwich)
library(pmdplyr)
library(stargazer)
library(DeclareDesign)
library(texreg)
library(metafor)
library(cjoint)
library(cregg)
library(ggtext)
library(xtable)
library
library(tidyverse)
library(ggpubr)

## set up working directories ##

onUser<-function(x){
  user<-Sys.info()["user"]
  onD<-grepl(x,user,ignore.case=TRUE)
  return(onD)
}
if(onUser("spb95")){
  setwd("~/Dropbox (Personal)/ypccc_parrish/covid_climate/replication")
  figfolder<-"~/Dropbox (Personal)/Apps/Overleaf/Covid Climate final BJPS/figures/"
  tabfolder<-"~/Dropbox (Personal)/Apps/Overleaf/Covid Climate final BJPS/tables/"
}
if(onUser("gabederoche")){
  setwd("~/Dropbox/covid_climate")
  figfolder<-"~/Dropbox/Apps/Overleaf/Covid Climate/figures/"
  tabfolder<-"~/Dropbox/Apps/Overleaf/Covid Climate/tables/"
}

cjtable<-function(d,caption,label){
  print.xtable(xtable(d,caption=caption),
               include.rownames=FALSE,
               file=paste0(tabfolder,label,".tex"),
               floating=TRUE,
               table.placement="h",
               label=paste0("tab:",label))
  }

## CONJOINT ANALYSIS #### 
### Canada Conjoint: main #### 
# ## Import CSV file from Qualtrics
# ## The file has two lines of string headers
ca <- read.qualtrics("donjoint_canada.csv",new.format=TRUE,
                     responses=c("Q260","Q266","Q269"),
                     covariates=c("Duration..in.seconds.","Q272_Page.Submit","Q248_Page.Submit","Q273_Page.Submit",
                                  # "Q243_1","Q243_2","Q243_3","Q243_4","Q243_5","Q243_6","Q243_7",
                                  "ideo_1"),
                     letter="G",
                     respondentID="ResponseId")

## re-order factor levels with first level as baseline
ca<-ca%>%
  mutate(ideology=case_when(ideo_1%in%c("0","1","2","3")~"left",
                            ideo_1%in%c("4","5","6")~"center",
                            ideo_1%in%c("7","8","9","10")~"right"),
         Q248_Page.Submit=as.numeric(Q248_Page.Submit),
         Q272_Page.Submit=as.numeric(Q272_Page.Submit),
         Q273_Page.Submit=as.numeric(Q273_Page.Submit),
         conjoint_timing=Q248_Page.Submit+Q272_Page.Submit+Q273_Page.Submit)
medtime<-median(ca$conjoint_timing,na.rm=TRUE)
ca<-ca%>%
  mutate(attn=ifelse(conjoint_timing>=medtime,1,0))
table(ca$ideology)
# table(ca$Q243_2)
# ca<-ca%>%
#   mutate(attn=ifelse(Q243_4==1&Q243_2==7,1,0)) ## attentive respondents are those who put "crime" in 1st and "unemployment" in 2nd position
table(ca$attn)

levels(ca$Assistance.to.corporations)
ca$Assistance.to.corporations<-factor(ca$Assistance.to.corporations,
                                      levels=c(levels(ca$Assistance.to.corporations)[3],
                                               levels(ca$Assistance.to.corporations)[1],
                                               levels(ca$Assistance.to.corporations)[2]),
                                      labels=c("No financial assistance for business",
                                               "Assistance for small businesses",
                                               "Assistance for corporations"))
levels(ca$Assistance.to.corporations)
levels(ca$Assistance.to.individuals)
ca$Assistance.to.individuals<-factor(ca$Assistance.to.individuals,
                                     levels=c(levels(ca$Assistance.to.individuals)[4],
                                              levels(ca$Assistance.to.individuals)[1],
                                              levels(ca$Assistance.to.individuals)[2],
                                              levels(ca$Assistance.to.individuals)[3]),
                                     labels=c("No financial assistance for individuals",
                                              "Universal basic income",
                                              "One-time cash payments for all",
                                              "Monthly cash payments if lost job"))
levels(ca$Assistance.to.individuals)
levels(ca$Climate.Action)
ca$Climate.Action<-factor(ca$Climate.Action,
                          levels=c(levels(ca$Climate.Action)[2],
                                   levels(ca$Climate.Action)[1],
                                   levels(ca$Climate.Action)[3],
                                   levels(ca$Climate.Action)[4]),
                          labels=c("No climate action",
                                   "Exclude fossil fuel companies",
                                   "Condition assistance on pollution reduction",
                                   "Retrain fossil fuel workers"))
levels(ca$Climate.Action)
levels(ca$Cost)
ca$Cost<-factor(ca$Cost,
                levels=c(levels(ca$Cost)[3],
                         levels(ca$Cost)[4],
                         levels(ca$Cost)[1],
                         levels(ca$Cost)[2]))
levels(ca$Cost)
levels(ca$Infrastructure.Investments)
ca$Infrastructure.Investments<-factor(ca$Infrastructure.Investments,
                                      levels=c(levels(ca$Infrastructure.Investments)[5],
                                               levels(ca$Infrastructure.Investments)[1],
                                               levels(ca$Infrastructure.Investments)[2],
                                               levels(ca$Infrastructure.Investments)[3],
                                               levels(ca$Infrastructure.Investments)[4]),
                                      labels=c("No infrastructure investments",
                                               "Clean energy","Clean transportation",
                                               "Fossil fuel infrastructure","Roads, bridges, and tunnels"))
levels(ca$Infrastructure.Investments)
ca<-ca%>%
  rename(`Business assistance`=Assistance.to.corporations,
         `Individual assistance`=Assistance.to.individuals,
         `Climate action`=Climate.Action,
         `Infrastructure`=Infrastructure.Investments)

amces.ca<-cj(ca,
             selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
             id=~respondent)

### Canada conjoint by left/right parties #### 
ca.l<-ca[ca$ideology=="left",]
id.l<-unique(ca.l$respondent)
id.l<-data.frame(respondent.old=id.l,respondent.new=1:length(id.l))
ca.l<-left_join(ca.l,id.l,by=c("respondent"="respondent.old"))
amces.l<-cj(ca.l,
            selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
            id=~respondent.new)

ca.c<-ca[ca$ideology=="center",]
id.c<-unique(ca.c$respondent)
id.c<-data.frame(respondent.old=id.c,respondent.new=1:length(id.c))
ca.c<-left_join(ca.c,id.c,by=c("respondent"="respondent.old"))
amces.c<-cj(ca.c,
            selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
            id=~respondent.new)

ca.r<-ca[ca$ideology=="right",]
id.r<-unique(ca.r$respondent)
id.r<-data.frame(respondent.old=id.r,respondent.new=1:length(id.r))
ca.r<-left_join(ca.r,id.r,by=c("respondent"="respondent.old"))
amces.r<-cj(ca.r,
            selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
            id=~respondent.new)

## 83% confidence intervals
amces.ca<-amces.ca%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$25 billion")==TRUE,"baseline","level"))
amces.r<-amces.r%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$25 billion")==TRUE,"baseline","level"))
amces.c<-amces.c%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$25 billion")==TRUE,"baseline","level"))

amces.l<-amces.l%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$25 billion")==TRUE,"baseline","level"))

labels<-amces.ca%>%
  filter(estimate==0)%>%
  mutate(type="label",level=feature)
amces.ca<-rbind(amces.ca,labels)
amces.r<-rbind(amces.r,labels)

## relevel for all 4 datasets. 
levels(amces.ca$level)
amces.ca$level<-factor(amces.ca$level,levels=c("$150 billion","$100 billion","$50 billion","$25 billion","Cost",
                                               "Universal basic income","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                               "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                               "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                               "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.c$level<-factor(amces.c$level,levels=c("$150 billion","$100 billion","$50 billion","$25 billion","Cost",
                                             "Universal basic income","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                             "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                             "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                             "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.r$level<-factor(amces.r$level,levels=c("$150 billion","$100 billion","$50 billion","$25 billion","Cost",
                                             "Universal basic income","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                             "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                             "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                             "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.l$level<-factor(amces.l$level,levels=c("$150 billion","$100 billion","$50 billion","$25 billion","Cost",
                                             "Universal basic income","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                             "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                             "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                             "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.r$ideology<-"Right"
amces.c$ideology<-"Center"
amces.l$ideology<-"Left"
amces.ca.ideo<-do.call(rbind,list(amces.r,amces.c,amces.l))

### US conjoint #### 
## add Q272_Page.Submit, Q248_Page.Submit, Q273_Page.Submit: these are the timers. 
us <- read.qualtrics("donjoint_us.csv",new.format=TRUE,
                     responses=c("Q260","Q266","Q269"),
                     covariates=c("Duration..in.seconds.","Q272_Page.Submit","Q248_Page.Submit","Q273_Page.Submit",
                                  # "Q243_1","Q243_2","Q243_3","Q243_4","Q243_5","Q243_6","Q243_7",
                                  "political_party","zip"),
                     letter="G",
                     respondentID="ResponseId")
us<-us%>%
  mutate(political_party=case_when(political_party==1~"Democrat",
                                   political_party==2~"Democrat",
                                   political_party==3~"Democrat",
                                   political_party==4~"Independent",
                                   political_party==5~"Republican",
                                   political_party==6~"Democrat",
                                   political_party==7~"Independent",
                                   political_party==8~"Republican",
                                   political_party==9~"Republican",
                                   political_party==10~"Republican"),
         Q248_Page.Submit=as.numeric(Q248_Page.Submit),
         Q272_Page.Submit=as.numeric(Q272_Page.Submit),
         Q273_Page.Submit=as.numeric(Q273_Page.Submit),
         conjoint_timing=Q248_Page.Submit+Q272_Page.Submit+Q273_Page.Submit)
medtime<-median(us$conjoint_timing,na.rm=TRUE)
us<-us%>%
  mutate(attn=ifelse(conjoint_timing>=medtime,1,0))
table(us$political_party)
table(us$attn)
## re-order factor levels with first level as baseline
levels(us$Assistance.to.corporations)
us$Assistance.to.corporations<-factor(us$Assistance.to.corporations,
                                      levels=c(levels(us$Assistance.to.corporations)[3],
                                               levels(us$Assistance.to.corporations)[2],
                                               levels(us$Assistance.to.corporations)[1]),
                                      labels=c("No financial assistance for business",
                                               "Assistance for small businesses",
                                               "Assistance for corporations"))
levels(us$Assistance.to.corporations)
levels(us$Assistance.to.individuals)
us$Assistance.to.individuals<-factor(us$Assistance.to.individuals,
                                     levels=c(levels(us$Assistance.to.individuals)[4],
                                              levels(us$Assistance.to.individuals)[2],
                                              levels(us$Assistance.to.individuals)[1],
                                              levels(us$Assistance.to.individuals)[3]),
                                     labels=c("No financial assistance for individuals",
                                              "Government healthcare for unemployed",
                                              "One-time cash payments for all",
                                              "Monthly cash payments if lost job"))
levels(us$Assistance.to.individuals)
levels(us$Climate.Action)
us$Climate.Action<-factor(us$Climate.Action,
                          levels=c(levels(us$Climate.Action)[2],
                                   levels(us$Climate.Action)[1],
                                   levels(us$Climate.Action)[3],
                                   levels(us$Climate.Action)[4]),
                          labels=c("No climate action",
                                   "Exclude fossil fuel companies",
                                   "Condition assistance on pollution reduction",
                                   "Retrain fossil fuel workers"))
levels(us$Climate.Action)
levels(us$Cost)
us$Cost<-factor(us$Cost,
                levels=c(levels(us$Cost)[4],
                         levels(us$Cost)[1],
                         levels(us$Cost)[2],
                         levels(us$Cost)[3]))
levels(us$Cost)
levels(us$Infrastructure.Investments)
us$Infrastructure.Investments<-factor(us$Infrastructure.Investments,
                                      levels=c(levels(us$Infrastructure.Investments)[5],
                                               levels(us$Infrastructure.Investments)[1],
                                               levels(us$Infrastructure.Investments)[2],
                                               levels(us$Infrastructure.Investments)[3],
                                               levels(us$Infrastructure.Investments)[4]),
                                      labels=c("No infrastructure investments",
                                               "Clean energy",
                                               "Clean transportation",
                                               "Fossil fuel infrastructure",
                                               "Roads, bridges, and tunnels"))
levels(us$Infrastructure.Investments)
us<-us%>%
  rename(`Business assistance`=Assistance.to.corporations,
         `Individual assistance`=Assistance.to.individuals,
         `Climate action`=Climate.Action,
         `Infrastructure`=Infrastructure.Investments)


amces.us<-cj(us,
             selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
             id=~respondent)

## AMCES by party 
us.R<-us[us$political_party=="Republican",]
id.R<-unique(us.R$respondent)
id.R<-data.frame(respondent.old=id.R,respondent.new=1:length(id.R))
us.R<-left_join(us.R,id.R,by=c("respondent"="respondent.old"))
amces.R<-cj(us.R,
            selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
            id=~respondent.new)

us.D<-us[us$political_party=="Democrat",]
id.D<-unique(us.D$respondent)
id.D<-data.frame(respondent.old=id.D,respondent.new=1:length(id.D))
us.D<-left_join(us.D,id.D,by=c("respondent"="respondent.old"))
amces.D<-cj(us.D,
            selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
            id=~respondent.new)

## make 83% confidence intervals for visual inspection of difference in attribute level effects 
qnorm((1-0.17/2)) ## multiply this by SE and add/ subtract from point estimates
amces.us<-amces.us%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion")==TRUE,"baseline","level"))
amces.R<-amces.R%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion")==TRUE,"baseline","level"))
amces.D<-amces.D%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion")==TRUE,"baseline","level"))


labels<-amces.us%>%
  filter(estimate==0)%>%
  mutate(type="label",level=feature)
amces.us<-rbind(amces.us,labels)
amces.R<-rbind(amces.R,labels)

amces.us$level<-factor(amces.us$level,levels=c("$3 trillion","$2 trillion","$1 trillion","$500 billion","Cost",
                                               "Government healthcare for unemployed","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                               "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                               "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                               "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
levels(amces.us$level)
amces.D$level<-factor(amces.D$level,levels=c("$3 trillion","$2 trillion","$1 trillion","$500 billion","Cost",
                                             "Government healthcare for unemployed","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                             "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                             "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                             "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.D$party<-"Democrat"
amces.R$level<-factor(amces.R$level,levels=c("$3 trillion","$2 trillion","$1 trillion","$500 billion","Cost",
                                             "Government healthcare for unemployed","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                             "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                             "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                             "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.R$party<-"Republican"

amces.pid<-rbind(amces.R,amces.D)

### conjoint plots #### 
bold.labels<-c("Infrastructure","Business assistance","Climate action","Individual assistance","Cost")
bold.labels2<-ifelse(levels(amces.us$level)%in%bold.labels,yes="bold",no="plain")
## US main plot 
plotus<-ggplot(amces.us,aes(x=level,y=estimate,shape=type,alpha=type#,color=feature
                            ))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.2)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.4)+
  scale_shape_manual(values=c(1,19,20),guide=FALSE)+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1,guide=FALSE))+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)",title="United States")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))
plotus
### Canada main plot
plotca<-ggplot(amces.ca,aes(x=level,y=estimate,shape=type,alpha=type#,color=feature
                            ))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.2)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.4)+
  scale_shape_manual(values=c(1,19,20),guide=FALSE)+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1,guide=FALSE))+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)",title="Canada")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))
plotca
## save these two plots together
ggsave(file=paste0(figfolder,"amces_us_ca.png"),width=8.5,height=7,units="in",
       ggarrange(plotus,plotca,nrow=1,ncol=2),dpi=600)

## US parties plot
amces.pid<-amces.pid%>%
  mutate(party=ifelse(type=="label"|type=="baseline","none",party))
plotus.pid<-ggplot(amces.pid,aes(x=level,y=estimate,shape=party,alpha=type))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.4)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.2)+
  scale_shape_manual(values=c("Republican"=17,"Democrat"=16,"none"=1),breaks=c("Republican","Democrat",NA))+ ## NA here removes 'none' from legend
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide="none")+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)")+
  ggtitle("US")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title=element_text(hjust=0.5))+
  guides(shape=guide_legend(nrow=2,byrow=TRUE))
## 
plotus.pid

## Canada left/ right ideology plot 
amces.ca.ideo<-amces.ca.ideo%>%
  mutate(ideology=ifelse(type=="label"|type=="baseline","none",ideology))
amces.ca.ideo$ideology<-factor(amces.ca.ideo$ideology,levels=c("Left","Center","Right","none"))
plotca.ideo<-ggplot(amces.ca.ideo,aes(x=level,y=estimate,shape=ideology,alpha=type))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.4)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.2)+
  scale_shape_manual(values=c("Left"=16,"Right"=17,"Center"=15,"none"=1),breaks=c("Left","Right","Center",NA))+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide="none")+
  # scale_color_manual(values=c("Left"="orange","Center"="red","Right"="blue"))+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.32,.32))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)",shape="")+
  ggtitle("Canada")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title=element_text(hjust=0.5))+
  guides(shape=guide_legend(nrow=2,byrow=TRUE))
plotca.ideo
## save ideology and partisan plot together 
ggsave(file=paste0(figfolder,"amces_ca_us_pid.png"),width=8.5,height=7,units="in",
       ggarrange(plotus.pid,plotca.ideo),dpi=600)

### conjoint tables ####
amces.us<-amces.us%>%
  select(feature,level,estimate,std.error,p)%>%
  rename(dimension=feature)%>%
  filter(!is.na(std.error))

cjtable(amces.us,"Full-sample results from conjoint experiment in the US","cjus_main")

amces.pid<-amces.pid%>%
  select(feature,level,estimate,std.error,p,party)%>%
  rename(dimension=feature)%>%
  filter(!is.na(std.error))
cjtable(amces.pid,"Results from conjoint experiment in the US by party","cjus_pid")

amces.ca<-amces.ca%>%
  select(feature,level,estimate,std.error,p)%>%
  rename(dimension=feature)%>%
  filter(!is.na(std.error))
cjtable(amces.ca,"Full-sample results from conjoint experiment in Canada","cjca_main")

amces.ca.ideo<-amces.ca.ideo%>%
  select(feature,level,estimate,std.error,p,ideology)%>%
  rename(dimension=feature)%>%
  filter(!is.na(std.error))
cjtable(amces.ca.ideo,"Results from conjoint experiment in Canada by ideology","cjca_ideo")
  
### Predicting support for hypothetical packages: Canada ####
ca2<-read_csv("conjoint_canada.csv")
selection<-names(ca2)[grep("G-",names(ca2))]### policy levels for each choice task --will only need the ones that have 7 characters--others are attribute labels
ca2<-ca2[3:nrow(ca2),]%>%
  select(ResponseId,all_of(selection),Q261_1,Q261_2,Q267_1,Q267_2,Q270_1,Q270_2,)
dimnames<-selection[which(nchar(selection)<=5)]

## record dimension names 
dimensions<-ca2%>%
  select(ResponseId,all_of(dimnames))%>%
  filter(!is.na(`G-1-1`))%>% ## these are incompletes 
  pivot_longer(cols=all_of(dimnames),values_to="dimname",names_to="dimid")%>%
  mutate(task=substr(dimid,start=3,stop=3),
         dimension=substr(dimid,start=5,stop=5))

levels<-selection[which(nchar(selection)>5)]
ca2<-ca2%>%
  select(-all_of(dimnames))%>%
  filter(!is.na(`G-1-1-1`))%>%
  pivot_longer(all_of(levels),values_to="level",names_to="levelid")%>%
  mutate(task=substr(levelid,start=3,stop=3),
         profile=substr(levelid,start=5,stop=5),
         dimension=substr(levelid,start=7,stop=7))
ca2<-ca2%>%
  left_join(dimensions,by=c("task","dimension","ResponseId"))

### manually add indicator for attribute to each level so that "nones" are distinguished in regression
ca2<-ca2%>%
  mutate(level=paste(dimname,level,sep="-"))

### set support column to correct value for each task x level combo
ca2<-ca2%>%
  mutate(support=case_when(task==1&profile==1~Q261_1,
                           task==1&profile==2~Q261_2,
                           task==2&profile==1~Q267_1,
                           task==2&profile==2~Q267_2,
                           task==3&profile==1~Q270_1,
                           task==3&profile==2~Q270_2))
unique(ca2$level)

### dummy out the levels ## 
ca2w<-ca2%>%
  mutate(v=1)%>%
  pivot_wider(id_cols=c(ResponseId,task,profile,dimension,support),names_from=level,values_from=v,values_fill=0)

ca2w<-ca2w[complete.cases(ca2w)==TRUE,]
ca2w$ResponseId<-factor(ca2w$ResponseId)
ca2w$support<-as.numeric(ca2w$support)
ca2w<-ca2w%>%
  select(-c("task","profile","dimension","Assistance to corporations-None",
            "Infrastructure Investments-None","Climate Action-None","Cost-$25 billion",
            "Assistance to individuals-None"))

## regression, omitting baseline categories for each level
l<-lm_robust(support~`Infrastructure Investments-Investments in clean energy, like wind and solar`+
               `Infrastructure Investments-Investments in clean transportation, like public transit and electric vehicle charging`+
               `Infrastructure Investments-Investments in facilities to generate electricity from fossil fuels (oil, gas, and coal)`+
               `Infrastructure Investments-Investments in roads, bridges, and tunnels`+
               `Assistance to corporations-Financial support to corporations if all workers keep their jobs`
             +`Assistance to corporations-Financial support to small businesses if all workers keep their jobs`
             +`Assistance to individuals-A guaranteed minimum annual income for every Canadian`
             +`Assistance to individuals-A one-time cash payment for every Canadian`
             +`Assistance to individuals-Monthly cash payments for every Canadian who has lost their job`
             +`Climate Action-Exclude oil, coal, and gas companies from receiving assistance`
             +`Climate Action-Require that companies reduce pollution in exchange for receiving financial assistance`
             +`Climate Action-Retrain workers in polluting industries, such as oil, coal, and gas`
             +`Cost-$100 billion`
             +`Cost-$150 billion`
             +`Cost-$50 billion`
             ,clusters=ResponseId,data=ca2w)

## predict values for 2 packages based on this regression 
p1.ca<-ca2w[1,]
names(p1.ca)

## set up data row for package which has lowest cost and includes no infrastructure spending or climate policy, 
## and includes assistance for corporations (b/c gets higher support than small business assistance) and monthly unemployment
## (for comparability with US, and b/c support is slightly higher-->conservative)
p1.ca[1,]<-list(NA,NA,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0)
predict(l,p1.ca,se.fit=TRUE) ## 61% support

## set up data row for package which holds individual and corporate assistance constant from above, 
## uses most popular climate policies (conditioning assistance on pollution reduction and 
## clean energy spending), and has highest cost 
p2.ca<-ca2w[1,]
p2.ca[1,]<-list(NA,NA,0,0,1,0,0,0,0,0,1,1,0,1,1,0,0)
predict(l,p2.ca,se.fit=TRUE) ## 75% support

### predicting support for hypothetical packages: US #### 
us2<-read_csv("conjoint_us.csv")
us2<-us2[3:nrow(us2),]%>%
  select(ResponseId,all_of(selection),Q261_1,Q261_2,Q267_1,Q267_2,Q270_1,Q270_2,)

## record dimension names
dimensions<-us2%>%
  select(ResponseId,all_of(dimnames))%>%
  filter(!is.na(`G-1-1`))%>% ## these are incompletes
  pivot_longer(cols=all_of(dimnames),values_to="dimname",names_to="dimid")%>%
  mutate(task=substr(dimid,start=3,stop=3),
         dimension=substr(dimid,start=5,stop=5))

us2<-us2%>%
  select(-all_of(dimnames))%>%
  filter(!is.na(`G-1-1-1`))%>%
  pivot_longer(all_of(levels),values_to="level",names_to="levelid")%>%
  mutate(task=substr(levelid,start=3,stop=3),
         profile=substr(levelid,start=5,stop=5),
         dimension=substr(levelid,start=7,stop=7))
us2<-us2%>%
  left_join(dimensions,by=c("task","dimension","ResponseId"))

### manually add indicator for attribute to each level so that "nones" are distinguished in regression
us2<-us2%>%
  mutate(level=paste(dimname,level,sep="-"))

### set support column to correct value for each task x level combo
us2<-us2%>%
  mutate(support=case_when(task==1&profile==1~Q261_1,
                           task==1&profile==2~Q261_2,
                           task==2&profile==1~Q267_1,
                           task==2&profile==2~Q267_2,
                           task==3&profile==1~Q270_1,
                           task==3&profile==2~Q270_2))
sort(unique(us2$level))

### dummy out the levels ##
us2w<-us2%>%
  mutate(v=1)%>%
  pivot_wider(id_cols=c(ResponseId,task,profile,dimension,support),names_from=level,values_from=v,values_fill=0)

us2w<-us2w[complete.cases(us2w)==TRUE,]
us2w$ResponseId<-factor(us2w$ResponseId)
us2w$support<-as.numeric(us2w$support)
us2w<-us2w%>%
  select(-c("task","profile","dimension","Assistance to corporations-None",
            "Infrastructure Investments-None","Climate Action-None","Cost-$500 billion",
            "Assistance to individuals-None"))

## regression, omitting baseline categories for each level
l2<-lm_robust(support~`Infrastructure Investments-Investments in clean energy, like wind and solar`+
                `Infrastructure Investments-Investments in clean transportation, like public transit and electric vehicle charging`+
                `Infrastructure Investments-Investments in facilities to generate electricity from fossil fuels (oil, gas, and coal)`+
                `Infrastructure Investments-Investments in roads, bridges, and tunnels`+
                `Assistance to corporations-Financial support to corporations if all workers keep their jobs`
              +`Assistance to corporations-Financial support to small businesses if all workers keep their jobs`
              +`Assistance to individuals-Government-sponsored healthcare for unemployed workers`
              +`Assistance to individuals-A one-time cash payment for every American`
              +`Assistance to individuals-Monthly cash payments for every American who has lost their job`
              +`Climate Action-Exclude oil, coal, and gas companies from receiving assistance`
              +`Climate Action-Require that companies reduce pollution in exchange for receiving financial assistance`
              +`Climate Action-Retrain workers in polluting industries, such as oil, coal, and gas`
              +`Cost-$1 trillion`
              +`Cost-$2 trillion`
              +`Cost-$3 trillion`
              ,clusters=ResponseId,data=us2w)

## next predict values for 2 packages based on this regression ##
## package with lowest cost and includes no infrastructure spending or climate policy, 
## and includes assistance for small businesses (real corp assistance option) and 
## monthly cash payments if lost job (real individual assistance--also indistinguishable from the other option)
p1.us<-us2w[1,]
p1.us[1,]<-list(NA,NA,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0)
predict(l2,p1.us,se.fit=TRUE) ## 66% support
p2.us<-us2w[1,]
## package which holds individual and corporate assistance constant, uses most popular climate policies (conditioning assistance on pollution reduction and
## clean energy spending), and has highest cost
p2.us[1,]<-list(NA,NA,0,1,0,1,1,0,0,0,0,1,0,0,0,1,0)
predict(l2,p2.us,se.fit=TRUE) ## 78% support

predtable<-xtable(do.call(rbind,list(
  unlist(predict(l,p1.ca,se.fit=TRUE)),
  unlist(predict(l,p2.ca,se.fit=TRUE)),
  unlist(predict(l2,p1.us,se.fit=TRUE)),
  unlist(predict(l2,p2.us,se.fit=TRUE)))),
  caption="Predicted support for hypothetical policy packages with and without climate policy elements. 
  The baseline packages (Rows 1 and 3) contain the most politically feasible COVID-19 response policies, 
  does not contain climate policy or infrastructure elements, and has the lowest cost. By comparison, 
  Rows 2 and 4 reflect support for packages with the most popular climate action and infrastructure elements 
  in each country and the highest cost. Predicted values reflect average support on a 0-100 scale. 
  Standard errors are clustered at the respondent level.", label="tab:predsupport")
colnames(predtable)<-c("Support","Std. Error")
rownames(predtable)<-c("Canada (baseline)","Canada (with climate policy)",
                       "US (baseline)","US (with climate policy)")

print.xtable(predtable, 
  file=paste0(tabfolder,"bundlecompare.tex"),
  floating=TRUE,
  table.placement="ht",include.rownames=TRUE
)

### follow-up conjoint with covid mitigation measures #### 
us2 <- read.qualtrics("conjoint2_us.csv",new.format=TRUE,
                      responses=c("con_choice_1","con_choice_2","con_choice_3"),
                      covariates=c("political_party","zip"),
                      letter="G",
                      respondentID="ResponseId")
## re-order factor levels with first level as baseline
levels(us2$Assistance.to.corporations)
us2$Assistance.to.corporations<-factor(us2$Assistance.to.corporations,
                                       levels=c(levels(us2$Assistance.to.corporations)[3],
                                                levels(us2$Assistance.to.corporations)[2],
                                                levels(us2$Assistance.to.corporations)[1]),
                                       labels=c("No financial assistance for business",
                                                "Assistance for small businesses",
                                                "Assistance for corporations"))
levels(us2$Assistance.to.corporations)
levels(us2$Assistance.to.individuals)
us2$Assistance.to.individuals<-factor(us2$Assistance.to.individuals,
                                      levels=c(levels(us2$Assistance.to.individuals)[4],
                                               levels(us2$Assistance.to.individuals)[2],
                                               levels(us2$Assistance.to.individuals)[1],
                                               levels(us2$Assistance.to.individuals)[3]),
                                      labels=c("No financial assistance for individuals",
                                               "Government healthcare for unemployed",
                                               "One-time cash payments for all",
                                               "Monthly cash payments if lost job"))
levels(us2$Assistance.to.individuals)
levels(us2$Climate.Action)
us2$Climate.Action<-factor(us2$Climate.Action,
                           levels=c(levels(us2$Climate.Action)[2],
                                    levels(us2$Climate.Action)[1],
                                    levels(us2$Climate.Action)[3],
                                    levels(us2$Climate.Action)[4]),
                           labels=c("No climate action",
                                    "Exclude fossil fuel companies",
                                    "Condition assistance on pollution reduction",
                                    "Retrain fossil fuel workers"))
levels(us2$Climate.Action)
levels(us2$Cost)
us2$Cost<-factor(us2$Cost,
                 levels=c(levels(us2$Cost)[4],
                          levels(us2$Cost)[1],
                          levels(us2$Cost)[2],
                          levels(us2$Cost)[3]))
levels(us2$Cost)
levels(us2$Infrastructure.Investments)
us2$Infrastructure.Investments<-factor(us2$Infrastructure.Investments,
                                       levels=c(levels(us2$Infrastructure.Investments)[5],
                                                levels(us2$Infrastructure.Investments)[1],
                                                levels(us2$Infrastructure.Investments)[2],
                                                levels(us2$Infrastructure.Investments)[3],
                                                levels(us2$Infrastructure.Investments)[4]),
                                       labels=c("No infrastructure investments",
                                                "Clean energy",
                                                "Clean transportation",
                                                "Fossil fuel infrastructure",
                                                "Roads, bridges, and tunnels"))
levels(us2$Infrastructure.Investments)
names(us2)
levels(us2$`COVID-19.management`)
us2$`COVID-19.management`<-factor(us2$`COVID-19.management`,
                                  levels=c(levels(us2$`COVID-19.management`)[4],
                                           levels(us2$`COVID-19.management`)[1],
                                           levels(us2$`COVID-19.management`)[2],
                                           levels(us2$`COVID-19.management`)[3]),
                                  labels=c("No COVID-19 management policy",
                                           "Funding for at-home tests",
                                           "National vaccine mandate",
                                           "National mask mandate"))
levels(us2$`COVID-19.management`)
us2<-us2%>%
  rename(`Business assistance`=Assistance.to.corporations,
         `Individual assistance`=Assistance.to.individuals,
         `Climate action`=Climate.Action,
         `Infrastructure`=Infrastructure.Investments,
         `COVID-19 management`=`COVID-19.management`)

amces.us2<-cj(us2,
              selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+`COVID-19 management`+Cost,
              id=~respondent)
amces.us2<-amces.us2%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion","No COVID-19 management policy")==TRUE,
                     "baseline","level"))
labels<-amces.us2%>%
  filter(estimate==0)%>%
  mutate(type="label",level=feature)
amces.us2<-rbind(amces.us2,labels)

unique(amces.us2$level)
amces.us2$level<-factor(amces.us2$level,levels=c("$3 trillion","$2 trillion","$1 trillion","$500 billion","Cost",
                                                 "Government healthcare for unemployed","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                                 "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                                 "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                                 "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure",
                                                 "Funding for at-home tests","National vaccine mandate","National mask mandate","No COVID-19 management policy","COVID-19 management"))
levels(amces.us2$level)

plotus2<-ggplot(amces.us2,aes(x=level,y=estimate,shape=type,alpha=type#,color=feature
                              ))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.2)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.4)+
  scale_shape_manual(values=c(1,19,20),guide=FALSE)+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1,guide=FALSE))+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        aspect.ratio=3,
        legend.position="none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

plotus2
ggsave(file=paste0(figfolder,"amces_us2.png"),width=8.5,height=7,units="in",
       plotus2,dpi=600)

## Conjoint robustness checks ####
### Check 1: results by attention #### 
#### US results by attention #### 
uslow<-us%>%
  filter(attn==0)
ushigh<-us%>%
  filter(attn==1)


idlow<-unique(uslow$respondent)
idlow<-data.frame(respondent.old=idlow,respondent.new=1:length(idlow))
uslow<-left_join(uslow,idlow,by=c("respondent"="respondent.old"))
amces.uslow<-cj(uslow,
            selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
            id=~respondent.new)
idhigh<-unique(ushigh$respondent)
idhigh<-data.frame(respondent.old=idhigh,respondent.new=1:length(idhigh))
ushigh<-left_join(ushigh,idhigh,by=c("respondent"="respondent.old"))
amces.ushigh<-cj(ushigh,
                selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
                id=~respondent.new)

amces.ushigh<-amces.ushigh%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion")==TRUE,"baseline","level"))
amces.uslow<-amces.uslow%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion")==TRUE,"baseline","level"))


labels<-amces.ushigh%>%
  filter(estimate==0)%>%
  mutate(type="label",level=feature)
amces.ushigh<-rbind(amces.ushigh,labels)
amces.uslow<-rbind(amces.uslow,labels)

amces.uslow$level<-factor(amces.uslow$level,levels=c("$3 trillion","$2 trillion","$1 trillion","$500 billion","Cost",
                                             "Government healthcare for unemployed","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                             "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                             "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                             "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.uslow$Attention<-"Low attention"
amces.ushigh$level<-factor(amces.ushigh$level,levels=c("$3 trillion","$2 trillion","$1 trillion","$500 billion","Cost",
                                             "Government healthcare for unemployed","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                             "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                             "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                             "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.ushigh$Attention<-"High attention"

amces.usattn<-rbind(amces.ushigh,amces.uslow)

#### Canada conjoint by attention ####
calow<-ca%>%
  filter(attn==0)
cahigh<-ca%>%
  filter(attn==1)


idlow<-unique(calow$respondent)
idlow<-data.frame(respondent.old=idlow,respondent.new=1:length(idlow))
calow<-left_join(calow,idlow,by=c("respondent"="respondent.old"))
amces.calow<-cj(calow,
                selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
                id=~respondent.new)
idhigh<-unique(cahigh$respondent)
idhigh<-data.frame(respondent.old=idhigh,respondent.new=1:length(idhigh))
cahigh<-left_join(cahigh,idhigh,by=c("respondent"="respondent.old"))
amces.cahigh<-cj(cahigh,
                 selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
                 id=~respondent.new)

amces.cahigh<-amces.cahigh%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion")==TRUE,"baseline","level"))
amces.calow<-amces.calow%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion")==TRUE,"baseline","level"))


labels<-amces.cahigh%>%
  filter(estimate==0)%>%
  mutate(type="label",level=feature)
amces.cahigh<-rbind(amces.cahigh,labels)
amces.calow<-rbind(amces.calow,labels)

amces.calow$level<-factor(amces.calow$level,levels=c("$150 billion","$100 billion","$50 billion","$25 billion","Cost",
                                                     "Universal basic income","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                                     "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                                     "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                                     "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.calow$Attention<-"Low attention"
amces.cahigh$level<-factor(amces.cahigh$level,levels=c("$150 billion","$100 billion","$50 billion","$25 billion","Cost",
                                                       "Universal basic income","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                                       "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                                       "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                                       "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
amces.cahigh$Attention<-"High attention"
amces.caattn<-rbind(amces.cahigh,amces.calow)

#### Plots by attention #### 
bold.labels<-c("Infrastructure","Business assistance","Climate action","Individual assistance","Cost")
bold.labels2<-ifelse(levels(amces.usattn$level)%in%bold.labels,yes="bold",no="plain")

## us plot by attention 
amces.usattn<-amces.usattn%>%
  mutate(Attention=ifelse(type=="label"|type=="baseline","none",Attention))
plotus.attn<-ggplot(amces.usattn,aes(x=level,y=estimate,shape=Attention,alpha=type))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.4)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.4)+
  scale_shape_manual(values=c("High attention"=17,"Low attention"=16,"none"=1),breaks=c("High attention","Low attention",NA))+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide="none")+
  # scale_color_manual(values=c("High attention"="orange","Low attention"="black"))+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)")+
  ggtitle("US")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title=element_text(hjust=0.5))
plotus.attn

bold.labels2<-ifelse(levels(amces.caattn$level)%in%bold.labels,yes="bold",no="plain")
amces.caattn<-amces.caattn%>%
  mutate(Attention=ifelse(type=="label"|type=="baseline","none",Attention))
plotca.attn<-ggplot(amces.caattn,aes(x=level,y=estimate,shape=Attention,alpha=type))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.4)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.4)+
  scale_shape_manual(values=c("High attention"=17,"Low attention"=16,"none"=1),breaks=c("High attention","Low attention",NA))+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide="none")+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)")+
  ggtitle("Canada")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title=element_text(hjust=0.5))
plotca.attn

ggsave(file=paste0(figfolder,"amces_byattn.png"),width=9.5,height=7,units="in",
       ggarrange(plotus.attn,plotca.attn,common.legend=TRUE),dpi=600)

### check 2: results if we exclude extremes #### 
## exclude packages with all "none's" and highest cost
table(ca$`Business assistance`)
table(ca$`Individual assistance`)
table(ca$`Climate action`)
table(ca$Cost)
table(ca$Infrastructure)

## set up variable showing packages with zero "none's" and lowest cost
## and with all 'none's' and highest cost

ca<-ca%>%
  mutate(e1=ifelse(`Business assistance`=="No financial assistance for business" & 
                     `Individual assistance`=="No financial assistance for individuals" & 
                     `Climate action`=="No climate action" & 
                     Cost =="$150 billion" & 
                     Infrastructure=="No infrastructure investments",1,0),
         e2=ifelse(`Business assistance`!="No financial assistance for business" & 
                     `Individual assistance`!="No financial assistance for individuals" & 
                     `Climate action`!="No climate action" & 
                     Cost =="$25 billion" & 
                     Infrastructure!="No infrastructure investments",1,0))
table(ca$e1) ## only 2 packages have all nones and highest cost
table(ca$e2) ## 526 packages have no nones and lowest cost--will run analysis with this. 

table(us$Cost)
us<-us%>%
  mutate(e1=ifelse(`Business assistance`=="No financial assistance for business" & 
                     `Individual assistance`=="No financial assistance for individuals" & 
                     `Climate action`=="No climate action" & 
                     Cost =="$3 trillion" & 
                     Infrastructure=="No infrastructure investments",1,0),
         e2=ifelse(`Business assistance`!="No financial assistance for business" & 
                     `Individual assistance`!="No financial assistance for individuals" & 
                     `Climate action`!="No climate action" & 
                     Cost =="$500 billion" & 
                     Infrastructure!="No infrastructure investments",1,0))
table(us$e1) ## 11 packages
table(us$e2) ## 728 packages

amces.car<-cj(ca[ca$e1==0&ca$e2==0,],
             selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
             id=~respondent)
amces.usr<-cj(us[us$e1==0&us$e2==0,],
              selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
              id=~respondent)
amces.usr<-amces.usr%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion")==TRUE,"baseline","level"))
labels<-amces.usr%>%
  filter(estimate==0)%>%
  mutate(type="label",level=feature)
amces.usr<-rbind(amces.usr,labels)

amces.usr$level<-factor(amces.usr$level,levels=c("$3 trillion","$2 trillion","$1 trillion","$500 billion","Cost",
                                               "Government healthcare for unemployed","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                               "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                               "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                               "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
levels(amces.usr$level)


amces.car<-amces.car%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$25 billion")==TRUE,"baseline","level"))

labels<-amces.car%>%
  filter(estimate==0)%>%
  mutate(type="label",level=feature)
amces.car<-rbind(amces.car,labels)
amces.car$level<-factor(amces.car$level,levels=c("$150 billion","$100 billion","$50 billion","$25 billion","Cost",
                                               "Universal basic income","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                               "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                               "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                               "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))


#### plots for conjoints excluding extremes #### 
bold.labels2<-ifelse(levels(amces.usr$level)%in%bold.labels,yes="bold",no="plain")

plotusr<-ggplot(amces.usr,aes(x=level,y=estimate,shape=type,alpha=type))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.2)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.4)+
  scale_shape_manual(values=c(1,19,20),guide=FALSE)+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1,guide=FALSE))+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)",title="United States")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))
plotusr

bold.labels2<-ifelse(levels(amces.car$level)%in%bold.labels,yes="bold",no="plain")
## US main plot 
plotcar<-ggplot(amces.car,aes(x=level,y=estimate,shape=type,alpha=type))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.2)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.4)+
  scale_shape_manual(values=c(1,19,20),guide=FALSE)+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1,guide=FALSE))+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)",title="Canada")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))
plotcar

ggsave(file=paste0(figfolder,"amces_woextremes.png"),width=9.5,height=7,units="in",
       ggarrange(plotusr,plotcar),dpi=600)

### check 3: split by round 
### check 3: split by round #### 
us$task<-factor(us$task)
diff_amcesus<-cj(us,
               selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
               id=~respondent,estimate="amce",by=~task)
ca$task<-factor(ca$task)
diff_amcesca<-cj(ca,
                 selected~`Business assistance`+`Individual assistance`+`Climate action`+Infrastructure+Cost,
                 id=~respondent,estimate="amce",by=~task)

#### plots for split by round #### 

diff_amcesus<-diff_amcesus%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$500 billion")==TRUE,"baseline","level"))

labels<-diff_amcesus%>%
  filter(estimate==0)%>%
  mutate(type="label",level=feature)
diff_amcesus<-rbind(diff_amcesus,labels)

diff_amcesus$level<-factor(diff_amcesus$level,levels=c("$3 trillion","$2 trillion","$1 trillion","$500 billion","Cost",
                                                     "Government healthcare for unemployed","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                                     "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                                     "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                                     "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
bold.labels2<-ifelse(levels(diff_amcesus$level)%in%bold.labels,yes="bold",no="plain")
diff_amcesus<-diff_amcesus%>%
  mutate(task=ifelse(type=="label"|type=="baseline","none",task))
plotus_rounds<-ggplot(diff_amcesus,aes(x=level,y=estimate,shape=task,alpha=type))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.2)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.4)+
  scale_shape_manual(values=c("1"=16,"2"=17,"3"=15,"none"=1),breaks=c("1","2","3",NA))+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide="none")+
  # scale_color_manual(values=c("1"="orange","2"="black","3"="green"))+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)",title="United States")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

plotus_rounds

diff_amcesca<-diff_amcesca%>%
  mutate(upper.83=estimate+1.37*std.error,
         lower.83=estimate-1.37*std.error,
         type=ifelse(level%in%c("No financial assistance for business","No financial assistance for individuals",
                                "No climate action","No infrastructure investments","$25 billion")==TRUE,"baseline","level"))

diff_amcesca<-rbind(diff_amcesca,labels)

diff_amcesca$level<-factor(diff_amcesca$level,levels=c("$150 billion","$100 billion","$50 billion","$25 billion","Cost",
                                                       "Universal basic income","One-time cash payments for all","Monthly cash payments if lost job","No financial assistance for individuals","Individual assistance",
                                                       "Exclude fossil fuel companies","Condition assistance on pollution reduction","Retrain fossil fuel workers","No climate action","Climate action",
                                                       "Assistance for corporations","Assistance for small businesses","No financial assistance for business","Business assistance",
                                                       "Roads, bridges, and tunnels","Fossil fuel infrastructure","Clean transportation","Clean energy","No infrastructure investments","Infrastructure"))
bold.labels2<-ifelse(levels(diff_amcesca$level)%in%bold.labels,yes="bold",no="plain")
diff_amcesca<-diff_amcesca%>%
  mutate(task=ifelse(type=="label"|type=="baseline","none",task))
plotca_rounds<-ggplot(diff_amcesca,aes(x=level,y=estimate,shape=task,alpha=type))+
  geom_errorbar(aes(ymin=lower.83,ymax=upper.83),position=position_dodge(0.75),size=0.4,width=0.2)+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(0.75),size=0.4)+
  scale_shape_manual(values=c("1"=16,"2"=17,"3"=15,"none"=1),breaks=c("1","2","3",NA))+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide="none")+
  scale_x_discrete(position="bottom")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.22,.22))+
  labs(x="",y="Effect on Pr(policy\npackage preferred)",title="Canada")+
  coord_flip()+
  theme_bw()+
  theme(axis.ticks.y=element_blank(),
        axis.text.y=element_text(face=bold.labels2),
        aspect.ratio=3,
        legend.position="bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))

plotca_rounds

ggsave(filename=paste0(figfolder,"amces_splitbyround.png"),width=9.5,height=7,units="in",
       ggarrange(plotus_rounds,plotca_rounds,common.legend=TRUE),dpi=600)
## PANEL DATA ANALYSIS ####
### DATA PREP: load data and recode ####

# Read in 5 wave panel data
wave5 <- read.csv("panel_5waves.csv")
        
# Add COVID prevalence data 
covdata <- read.csv("covprev.csv")
#Merge to panel data by respondent ID
data <- left_join(wave5, covdata, by=c("responseid", "wave"))

remove(covdata)
remove(wave5)

#Recode time index
data$wave <- recode(data$wave, "wave1" = 1, "wave2" = 2, "wave3" = 3, "wave4" = 4, "wave5" = 5)

#Keep only waves 4 and 5
data <- subset(data, data$wave==4 | data$wave==5)

#covid prevalence = 0.000001 for wave 4 (so that we can log this)
data$covprev[data$wave == 4] <- 0

# COVID prevalence as a percentage rather than per 100,000
data$covprev <- data$covprev / 1000

#### Issue importance variables

#Climate
data$issue_climate <- coalesce(data$p2_5, data$P2_5r5)
data$issue_climate <- recode(data$issue_climate, "Extremely important" = 5, "Extremely important1" = 5, "Very important" = 4, "Very important2" = 4, "Somewhat important" = 3, "Somewhat important3" = 3, "Not too important" = 2, "Not too important4" = 2, "Not at all important" = 1, "Not at all important5" = 1)

#Jobs & Economy
data$issue_econ <- coalesce(data$p2_1, data$P2_5r1)
data$issue_econ <- recode(data$issue_econ, "Extremely important" = 5, "Extremely important1" = 5, "Very important" = 4, "Very important2" = 4, "Somewhat important" = 3, "Somewhat important3" = 3, "Not too important" = 2, "Not too important4" = 2, "Not at all important" = 1, "Not at all important5" = 1)

#Health care
data$issue_health <- coalesce(data$p2_2, data$P2_5r2)
data$issue_health <- recode(data$issue_health, "Extremely important" = 5, "Extremely important1" = 5, "Very important" = 4, "Very important2" = 4, "Somewhat important" = 3, "Somewhat important3" = 3, "Not too important" = 2, "Not too important4" = 2, "Not at all important" = 1, "Not at all important5" = 1)

#Education
data$issue_educ <- coalesce(data$p2_3, data$P2_5r3)
data$issue_educ <- recode(data$issue_educ, "Extremely important" = 5, "Extremely important1" = 5, "Very important" = 4, "Very important2" = 4, "Somewhat important" = 3, "Somewhat important3" = 3, "Not too important" = 2, "Not too important4" = 2, "Not at all important" = 1, "Not at all important5" = 1)

#Taxes and govt spending
data$issue_taxspend <- coalesce(data$p2_4, data$P2_5r4)
data$issue_taxspend <- recode(data$issue_taxspend, "Extremely important" = 5, "Extremely important1" = 5, "Very important" = 4, "Very important2" = 4, "Somewhat important" = 3, "Somewhat important3" = 3, "Not too important" = 2, "Not too important4" = 2, "Not at all important" = 1, "Not at all important5" = 1)

#Affordability
data$issue_afford <- coalesce(data$p2_6, data$P2_5r6)
data$issue_afford <- recode(data$issue_afford, "Extremely important" = 5, "Extremely important1" = 5, "Very important" = 4, "Very important2" = 4, "Somewhat important" = 3, "Somewhat important3" = 3, "Not too important" = 2, "Not too important4" = 2, "Not at all important" = 1, "Not at all important5" = 1)

#Immigration
data$issue_immig <- coalesce(data$p2_7, data$P2_5r7)
data$issue_immig <- recode(data$issue_immig, "Extremely important" = 5, "Extremely important1" = 5, "Very important" = 4, "Very important2" = 4, "Somewhat important" = 3, "Somewhat important3" = 3, "Not too important" = 2, "Not too important4" = 2, "Not at all important" = 1, "Not at all important5" = 1)

### Carbon pricing DVs

#Carbon Pricing Support
data$cp_support <- coalesce(data$po2, data$PO2_5)
data$cp_support <- recode(data$cp_support, "Strongly oppose" = 1, "Somewhat oppose" = 2, "Not sure" = 3, "Somewhat support" = 4, "Strongly support" = 5)

# Fair for average Canadians 
data$po10_1 <- recode(data$po10_1, "Strongly disagree1"=1, "2"=2,"Neither agree nor disagree3"=3, "4"=4, "Strongly agree5"=5)
data$cp_fairavg <- coalesce(data$po10_1, data$PO10_5r1)

# Negative impact on economy
data$cp_negimpact <- NA
data$cp_negimpact[data$wave==4] <- data$po10_5[data$wave==4]
data$cp_negimpact[data$wave==5] <- data$PO10_5r5[data$wave==5]
table(data$cp_negimpact, useNA="always")
data$cp_negimpact <- recode(data$cp_negimpact, "Strongly disagree1"=1, "1"=1, "2"=2, "Neither agree nor disagree3"=3, "3"=3, "4"=4, "Strongly agree5"=5, "5"=5)
class(data$cp_negimpact)

### DATA PREP FOR TWO-WAY FE ANALYSIS ####

# Keep only the columns with IV & DV 
twfedata <- data %>% 
  select(responseid, wave, 
         issue_climate, issue_econ, issue_health, issue_educ, issue_taxspend, issue_afford, issue_immig,
         cp_support, cp_fairavg, cp_negimpact, 
         covprev)

####Make the panel balanced
#Remove units with no covprev data
incompleteids <- twfedata$responseid[which(is.na(twfedata$covprev))]
twfedata <- twfedata[-which(twfedata$responseid %in% incompleteids),]
#Remove units with no wave 5 measures
twfedata <- make.pbalanced(twfedata, balance.type="shared.individuals", index=c("responseid","wave"))

remove(incompleteids)

### TWFE ANALYSIS ####
### TWFE: ISSUE IMPORTANCE ###########
twfedata<-mutate(twfedata,lncovprev=log(covprev+1)) ## add one here so that we can log the variable
#Issue - Climate
issue_climate_mod <- plm(formula=issue_climate ~ lncovprev, 
                         data=twfedata, model="within", effect="twoways",
                         index=c("responseid", "wave"))
coeftest(issue_climate_mod, vcov=vcovHC(issue_climate_mod, type="HC1", cluster="group"))

#Issue - Jobs and Economy
issue_econ_mod <- plm(formula=issue_econ ~ lncovprev, 
                      data=twfedata, model="within", effect="twoways",
                      index=c("responseid", "wave"))
coeftest(issue_econ_mod, vcov=vcovHC(issue_econ_mod, type="HC1", cluster="group"))
summary(issue_climate_mod) ## coef = .06 (166), p=.7

#Issue - Health Care
issue_health_mod <- plm(formula=issue_health ~ lncovprev, 
                        data=twfedata, model="within", effect="twoways",
                        index=c("responseid", "wave"))
coeftest(issue_health_mod, vcov=vcovHC(issue_health_mod, type="HC1", cluster="group"))
summary(issue_health_mod) ## coef = .03 (.13), p=.8

#Issue - Education 
issue_educ_mod <- plm(formula=issue_educ ~ lncovprev, 
                      data=twfedata, model="within", effect="twoways",
                      index=c("responseid", "wave"))
coeftest(issue_educ_mod, vcov=vcovHC(issue_educ_mod, type="HC1", cluster="group"))
summary(issue_educ_mod) ## coef = .2 (.18), p=.27

#Issue - Taxes and Govt Spending 
issue_taxspend_mod <- plm(formula=issue_taxspend ~ lncovprev, 
                          data=twfedata, model="within", effect="twoways",
                          index=c("responseid", "wave"))
coeftest(issue_taxspend_mod, vcov=vcovHC(issue_taxspend_mod, type="HC1", cluster="group"))
summary(issue_taxspend_mod) ## coef = -0.12 (.16), p=.47

#Issue - Affordability 
issue_afford_mod <- plm(formula=issue_afford ~ lncovprev, 
                        data=twfedata, model="within", effect="twoways",
                        index=c("responseid", "wave"))
coeftest(issue_afford_mod, vcov=vcovHC(issue_afford_mod, type="HC1", cluster="group"))
summary(issue_afford_mod) ## coef = -0.14 (0.15), p=.35

#Issue - Immigration 
issue_immig_mod <- plm(formula=issue_immig ~ lncovprev, 
                       data=twfedata, model="within", effect="twoways",
                       index=c("responseid", "wave"))
coeftest(issue_immig_mod, vcov=vcovHC(issue_immig_mod, type="HC1", cluster="group"))
summary(issue_immig_mod) ## coef = .04 (0..22), p=.87

### TWFE ANALYSIS - CARBON PRICING DVs ###########

##Support for Carbon Pricing 
#Balance the panel for this DV (only one where there is missing data) 
incids_cpsupport <- twfedata$responseid[which(is.na(twfedata$cp_support))]
twfedata_cpsupport <- twfedata[-which(twfedata$responseid %in% incids_cpsupport),]
#Model with these data only
cp_support_mod <- plm(formula= cp_support ~ lncovprev, 
                      data=twfedata_cpsupport, model="within", effect="twoways",
                      index=c("responseid", "wave"))
coeftest(cp_support_mod, vcov=vcovHC(cp_support_mod, type="HC1", cluster="group"))
summary(cp_support_mod) ## coef = -0.13 (0.22), p=.58

remove(incids_cpsupport)

##Fair for average Cdns (back to full dataset)
cp_fairavg_mod <- plm(formula= cp_fairavg ~ lncovprev, 
                      data=twfedata, model="within", effect="twoways",
                      index=c("responseid", "wave"))
coeftest(cp_fairavg_mod, vcov=vcovHC(cp_fairavg_mod, type="HC1", cluster="group"))
summary(cp_fairavg_mod) ## coef = .21 (.22), p=.35

##Negative impact on economy (full dataset)
cp_negimpact_mod <- plm(formula= cp_negimpact ~ lncovprev, 
                        data=twfedata, model="within", effect="twoways",
                        index=c("responseid", "wave"))
coeftest(cp_negimpact_mod, vcov=vcovHC(cp_negimpact_mod, type="HC1", cluster="group"))
summary(cp_negimpact_mod) ## coef = .03 (.24), p=.89

### Export Table 1 in Paper ####

#Get robust standard errors for regression table 
issue_climate_cov <- vcovHC(issue_climate_mod, type="HC1", cluster="group")
issue_climate_se <- sqrt(diag(issue_climate_cov))
cp_support_cov <- vcovHC(cp_support_mod, type="HC1", cluster="group")
cp_support_se <- sqrt(diag(cp_support_cov))

#Produce table
stargazer(issue_climate_mod, cp_support_mod,
          title = "Results of two-way fixed effects models, showing no statistically significant effect of COVID-19 prevalence at the local level on (1) the importance Canadians assign to climate change and (2) support for carbon pricing. Local COVID-19 prevalence is measured as the natural log of the percent of the local population that had contracted COVID-19 (cumulative through May 18, 2020). Issue importance is measured on a five-point scale (1:5). Support for carbon pricing is measured on a five-point scale (1:5)",
          label = "tab:clim_issueimp_nonstandardized",
          covariate.labels = c("COVID-19 prevalence (log)"),
          dep.var.labels = c("Importance of clim. chg.", "Support for carbon pricing"),
          omit = c("Constant"),
          model.names = FALSE,
          add.lines = list(c("Respondent fixed effects", "Yes", "Yes"),
                           c("Wave fixed effects", "Yes", "Yes"),
                           c("Waves", "2", "2"),
                           c("Observations", "834", "505")),
          se = list(issue_climate_se, cp_support_se),
          omit.stat = c("all"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = "Robust standard errors clustered at the respondent level",
          digits = 3,
          initial.zero = FALSE,
          header = FALSE,table.placement="!h",
          out = paste0(tabfolder,"tab1_twfe_results.tex"))


### SCALED TWFE ANALYSIS ######

# Scale all issue and climate policy variables
twfedata <- twfedata %>%
  mutate(issue_climate_sc = scale(issue_climate),
         issue_econ_sc = scale(issue_econ),
         issue_health_sc = scale(issue_health),
         issue_educ_sc = scale(issue_econ),
         issue_taxspend_sc = scale(issue_taxspend),
         issue_afford_sc = scale(issue_afford),
         issue_immig_sc = scale(issue_immig),
         cp_support_sc = scale(cp_support),
         cp_fairavg_sc = scale(cp_fairavg),
         cp_negimpact_sc = scale(cp_negimpact))

#Issue: Climate
issue_climate_mod_sc <- plm(formula=issue_climate_sc ~ lncovprev,
                            data=twfedata, model="within", effect="twoways",
                            index=c("responseid", "wave"))
coeftest(issue_climate_mod_sc, vcov=vcovHC(issue_climate_mod_sc, type="HC1", cluster="group"))

#Issue: Jobs & Economy
issue_econ_mod_sc <- plm(formula=issue_econ_sc ~ lncovprev,
                            data=twfedata, model="within", effect="twoways",
                            index=c("responseid", "wave"))
coeftest(issue_econ_mod_sc, vcov=vcovHC(issue_econ_mod_sc, type="HC1", cluster="group"))

#Issue: Health care
issue_health_mod_sc <- plm(formula=issue_health_sc ~ lncovprev,
                         data=twfedata, model="within", effect="twoways",
                         index=c("responseid", "wave"))
coeftest(issue_health_mod_sc, vcov=vcovHC(issue_health_mod_sc, type="HC1", cluster="group"))

#Issue: Education
issue_educ_mod_sc <- plm(formula=issue_educ_sc ~ lncovprev,
                           data=twfedata, model="within", effect="twoways",
                           index=c("responseid", "wave"))
coeftest(issue_educ_mod_sc, vcov=vcovHC(issue_educ_mod_sc, type="HC1", cluster="group"))

#Issue: Taxes and government spending
issue_taxspend_mod_sc <- plm(formula=issue_taxspend_sc ~ lncovprev,
                         data=twfedata, model="within", effect="twoways",
                         index=c("responseid", "wave"))
coeftest(issue_taxspend_mod_sc, vcov=vcovHC(issue_taxspend_mod_sc, type="HC1", cluster="group"))

#Issue: Affordability
issue_afford_mod_sc <- plm(formula=issue_afford_sc ~ lncovprev,
                             data=twfedata, model="within", effect="twoways",
                             index=c("responseid", "wave"))
coeftest(issue_afford_mod_sc, vcov=vcovHC(issue_afford_mod_sc, type="HC1", cluster="group"))

#Issue: Immigration
issue_immig_mod_sc <- plm(formula=issue_immig_sc ~ lncovprev,
                           data=twfedata, model="within", effect="twoways",
                           index=c("responseid", "wave"))
coeftest(issue_immig_mod_sc, vcov=vcovHC(issue_immig_mod_sc, type="HC1", cluster="group"))

### Climate Policy DVs

##Support for Carbon Pricing 
#Balance the panel for this DV (only one where there is missing data) 
incids_cpsupport <- twfedata$responseid[which(is.na(twfedata$cp_support_sc))]
twfedata_cpsupport <- twfedata[-which(twfedata$responseid %in% incids_cpsupport),]
#Model with these data only
cp_support_mod_sc <- plm(formula= cp_support_sc ~ lncovprev, 
                      data=twfedata_cpsupport, model="within", effect="twoways",
                      index=c("responseid", "wave"))
coeftest(cp_support_mod_sc, vcov=vcovHC(cp_support_mod_sc, type="HC1", cluster="group"))
remove(incids_cpsupport)

##Fair for avg. Cdn. 
cp_fairavg_mod_sc <- plm(formula= cp_fairavg_sc ~ lncovprev, 
                      data=twfedata, model="within", effect="twoways",
                      index=c("responseid", "wave"))
coeftest(cp_fairavg_mod_sc, vcov=vcovHC(cp_fairavg_mod_sc, type="HC1", cluster="group"))

##Fair for avg. Cdn. 
cp_negimpact_mod_sc <- plm(formula= cp_negimpact_sc ~ lncovprev, 
                         data=twfedata, model="within", effect="twoways",
                         index=c("responseid", "wave"))
coeftest(cp_negimpact_mod_sc, vcov=vcovHC(cp_negimpact_mod_sc, type="HC1", cluster="group"))

### Export Table 1 (SCALED) for Paper #### 

#Get robust standard errors for regression table 
issue_climate_cov_sc <- vcovHC(issue_climate_mod_sc, type="HC1", cluster="group")
issue_climate_se_sc <- sqrt(diag(issue_climate_cov_sc))
cp_support_cov_sc <- vcovHC(cp_support_mod_sc, type="HC1", cluster="group")
cp_support_se_sc <- sqrt(diag(cp_support_cov_sc))

#Produce table
stargazer(issue_climate_mod_sc, cp_support_mod_sc,
          title = "Results of two-way fixed effects models, showing no statistically significant effect of COVID-19 prevalence at the local level on (1) the importance Canadians assign to climate change and (2) support for carbon pricing. Local COVID-19 prevalence is measured as the natural log of the percent of the local population that had contracted COVID-19 (cumulative through May 18, 2020). Issue importance is measured on a five-point scale (1:5), but has been rescaled to have mean of zero and standard-deviation of one. Support for carbon pricing is measured on a five-point scale (1:5), but has also been rescaled to have mean of zero and standard-deviation of one. The table reports standardized coefficients, such that the effects should be interpreted in standard-deviation units.",
          label = "tab:clim_issueimp",
          covariate.labels = c("COVID-19 prevalence (log)"),
          dep.var.labels = c("Importance of clim. chg. (st.dev.)", "Support for carbon pricing (st.dev.)"),
          omit = c("Constant"),
          model.names = FALSE,
          add.lines = list(c("Respondent fixed effects", "Yes", "Yes"),
                           c("Wave fixed effects", "Yes", "Yes"),
                           c("Waves", "2", "2"),
                           c("Observations", "834", "505")),
          se = list(issue_climate_se_sc, cp_support_se_sc),
          omit.stat = c("all"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = "Robust standard errors clustered at the respondent level",
          digits = 3,
          initial.zero = FALSE,
          header = FALSE,
          out = paste0(tabfolder,"tab1_twfe_results_standardized.tex"))

### MDEs for panel analysis #### 
str(model.matrix(issue_climate_mod_sc)) ## N=1,668
summary(cp_support_mod_sc) ## N=1,010
# responseid<-seq(1,834,1)
# Z=simple_ra(834,prob=0.5,conditions=c("control","treatment"))
# y0<-rnorm(834,mean=0,sd=1)
# effect_size=0.5
# y1<-y0+effect_size*(Z=="treatment")

## MDE function: 
mdefunction<-function(sims,N,effectsizes=seq(0.01,0.35,0.01)){
  pval<-list()
for(i in 1:length(effectsizes)){
  print(i)
  effect_size<-effectsizes[i]
  pval[[i]]<-matrix(NA,nrow=1,ncol=sims)
for(j in 1:sims){
d<-data.frame(responseid=factor(rep(seq(1,N/2,1),2)),
  Z=simple_ra(N,prob=0.5,conditions=c("control","treatment")),
  y0=rnorm(N,mean=0,sd=1))%>%
    mutate(y1=y0+effect_size*(Z=="treatment")
  )
d$wave<-c(rep(4,N/2),rep(5,N/2))
d$wave<-factor(d$wave)
lm1<-plm(formula= y1 ~ Z,data=d, model="within", effect="twoways",
                              index=c("responseid","wave")) ### not sure how to handle the wave FEs here. 
pval[[i]][,j]<-coeftest(lm1, vcov=vcovHC(lm1, type="HC1", cluster="group"))[4]<0.05
}
}
return(pval)
}

## MDE for issue priority
pwr<-mdefunction(1000,1668)
pwr<-do.call(rbind,pwr)
sizes<-seq(0.01,0.35,0.01)
MDE<-cbind(sizes,apply(pwr,MARGIN=1,FUN=mean))
MDE
MDE<-data.frame(MDE)
names(MDE)<-c("effect","power")
write.csv(MDE,file="figures/pwr_panel_issuesupport.csv",row.names = FALSE)

## MDE for carbon pricing support
pwr2<-mdefunction(1000,1010,effectsizes=seq(0.01,0.35,0.01))
pwr2<-do.call(rbind,pwr2)
sizes<-seq(0.01,0.35,0.01)
MDE2<-cbind(sizes,apply(pwr2,MARGIN=1,FUN=mean))
MDE2<-data.frame(MDE2)
names(MDE2)<-c("effect","power")
write.csv(MDE2,file="figures/pwr_panel_cpsupport.csv",row.names = FALSE)

### figures for MDEs from panel analysis #### 
## load saved power analysis data 
MDE<-read_csv("figures/pwr_panel_issuesupport.csv")
MDE2<-read_csv("figures/pwr_panel_cpsupport.csv")


ggsave(filename=paste0(figfolder,"powercurve_issuepriority.png"),width=6,height=4,dpi=600,
       ggplot(MDE,aes(x=effect,y=power))+
         geom_point()+
         geom_smooth(se=FALSE,color="black",lwd=0.5)+
         geom_hline(yintercept=0.8,lty="dashed",color="darkgray")+
         labs(x="Effect size (std.dev)",y="Power")+
         theme_bw()
)

ggsave(filename=paste0(figfolder,"powercurve_carbonprice.png"),width=6,height=4,dpi=600,
ggplot(MDE2,aes(x=effect,y=power))+
  geom_point()+
  geom_smooth(se=FALSE,color="black",lwd=0.5)+
  geom_hline(yintercept=0.8,lty="dashed",color="darkgray")+
  labs(x="Effect size (std.dev)",y="Power")+
  theme_bw()
)


## SURVEY EXPERIMENTS #### 
### Canada 2021 ##########
#Read in Canada 2021 survey data
cansvy <- read.csv("conjoint_canada.csv", na.strings = c("", " "))

#Drop observations who didn't give consent
cansvy <- subset(cansvy, cansvy$Intro=="Yes")
#Drop observations who failed attention/manip check 
cansvy <- subset(cansvy, cansvy$expcheck=="Exponential growth" | cansvy$exprecheck=="Exponential growth")

#Recode DVs 
cansvy$climateworry <- recode(cansvy$climateworry, "Not at all worried"=0, "Not very worried"=1, "Somewhat worried"=2, "Very worried"=3)
cansvy$climatepriority <- recode(cansvy$climatepriority, "Low"=0, "Medium"=1, "High"=2, "Very high"=3)

#Indicators for treatment groups
cansvy$climexp_cont <- ifelse(cansvy$climateexperiment=="A", 1, 0)
cansvy$climexp_t_txt <- ifelse(cansvy$climateexperiment=="B", 1,
                               ifelse(cansvy$climateexperiment=="A", 0, NA))
cansvy$climexp_t_img <- ifelse(cansvy$climateexperiment=="C", 1, 
                               ifelse(cansvy$climateexperiment=="A", 0, NA))
cansvy$climexp_trt <- ifelse(cansvy$climateexperiment=="A", 0, 1)

#### ATEs

## Climate worry:

#Text-only treatment
ate_worry_can_txt <- lm(climateworry ~ climexp_t_txt, data=cansvy)
coeftest(ate_worry_can_txt, vcov=vcovHC(ate_worry_can_txt, type="HC1"))
#Image treatment effect
ate_worry_can_img <- lm(climateworry ~ climexp_t_img, data=cansvy)
coeftest(ate_worry_can_img, vcov=vcovHC(ate_worry_can_img, type="HC1"))  
#Any treatment
ate_worry_can_trt <- lm(climateworry ~ climexp_trt, data=cansvy)
coeftest(ate_worry_can_trt, vcov=vcovHC(ate_worry_can_trt, type="HC1"))

## Climate priority:

#Text-only treatment
ate_priority_can_txt <- lm(climatepriority ~ climexp_t_txt, data=cansvy)
coeftest(ate_priority_can_txt, vcov=vcovHC(ate_priority_can_txt, type="HC1"))
#Image treatment
ate_priority_can_img <- lm(climatepriority ~ climexp_t_img, data=cansvy)
coeftest(ate_priority_can_img, vcov=vcovHC(ate_priority_can_img, type="HC1"))
#Any treatment
ate_priority_can_trt <- lm(climatepriority ~ climexp_trt, data=cansvy)
coeftest(ate_priority_can_trt, vcov=vcovHC(ate_priority_can_trt, type="HC1"))

### USA 2021 ##########

# Read in 2021 US survey data
usasvy <- read.csv("conjoint_us.csv")
#Drop observations that didn't give consent
usasvy <- subset(usasvy, usasvy$Intro=="Yes")
#Drop observations who failed first pair of attention checks
usasvy <- subset(usasvy, usasvy$expcheck=="Exponential growth" | usasvy$exprecheck=="Exponential growth")

#Recode DVs 
usasvy$climateworry <- recode(usasvy$climateworry, "Not at all worried"=0, "Not very worried"=1, "Somewhat worried"=2, "Very worried"=3)
usasvy$climatepriority <- recode(usasvy$climatepriority, "Low"=0, "Medium"=1, "High"=2, "Very high"=3)

#Indicators for treatment groups
usasvy$climexp_cont <- ifelse(usasvy$climateexperiment=="A", 1, 0)
usasvy$climexp_t_txt <- ifelse(usasvy$climateexperiment=="B", 1,
                               ifelse(usasvy$climateexperiment=="A", 0, NA))
usasvy$climexp_t_img <- ifelse(usasvy$climateexperiment=="C", 1, 
                               ifelse(usasvy$climateexperiment=="A", 0, NA))
usasvy$climexp_trt <- ifelse(usasvy$climateexperiment=="A", 0, 1)

#### ATEs

## Climate worry:

#Text-only treatment
ate_worry_usa_txt <- lm(climateworry ~ climexp_t_txt, data=usasvy)
coeftest(ate_worry_usa_txt, vcov=vcovHC(ate_worry_usa_txt, type="HC1"))
#Image treatment effect
ate_worry_usa_img <- lm(climateworry ~ climexp_t_img, data=usasvy)
coeftest(ate_worry_usa_img, vcov=vcovHC(ate_worry_usa_img, type="HC1"))  
#Any treatment
ate_worry_usa_trt <- lm(climateworry ~ climexp_trt, data=usasvy)
coeftest(ate_worry_usa_trt, vcov=vcovHC(ate_worry_usa_trt, type="HC1"))

## Climate Priority:

#Text-only treatment
ate_priority_usa_txt <- lm(climatepriority ~ climexp_t_txt, data=usasvy)
coeftest(ate_priority_usa_txt, vcov=vcovHC(ate_priority_usa_txt, type="HC1"))
#Image treatment
ate_priority_usa_img <- lm(climatepriority ~ climexp_t_img, data=usasvy)
coeftest(ate_priority_usa_img, vcov=vcovHC(ate_priority_usa_img, type="HC1"))
#Any treatment
ate_priority_usa_trt <- lm(climatepriority ~ climexp_trt, data=usasvy)
coeftest(ate_priority_usa_trt, vcov=vcovHC(ate_priority_usa_trt, type="HC1"))


### Canada 2020 ##########

#Read in Canada 2020 survey data
can2020 <- data

# Code treatment groups 
can2020$climexp_cont <- ifelse(can2020$GROUP_EXP=="Group 1", 1, 0)   
can2020$climexp_t_txt <- ifelse(can2020$GROUP_EXP=="Group 2", 1, 
                                ifelse(can2020$GROUP_EXP=="Group 1", 0, NA))
can2020$climexp_t_img <- ifelse(can2020$GROUP_EXP=="Group 3", 1, 
                                ifelse(can2020$GROUP_EXP=="Group 1", 0, NA))
can2020$climexp_trt <- ifelse(can2020$GROUP_EXP=="Group 2" | can2020$GROUP_EXP=="Group 3", 1, 0)

# DV - only have climate worry (not climate priority in the 2020 data)
can2020$climateworry <- can2020$A3_5
can2020$climateworry <- na_if(can2020$climateworry, "Not sure")
can2020$climateworry <- recode(can2020$climateworry, "Not at all worried"=0, "Not too worried"=1, "Somewhat worried"=2, "Very worried"=3)

# 2020 Canada ATEs
ate_worry_can2020_txt <- lm(climateworry ~ climexp_t_txt, data=can2020)
coeftest(ate_worry_can2020_txt, vcov=vcovHC(ate_worry_can2020_txt, type="HC1"))
ate_worry_can2020_img <- lm(climateworry ~ climexp_t_img, data=can2020)
coeftest(ate_worry_can2020_img, vcov=vcovHC(ate_worry_can2020_img, type="HC1"))
ate_worry_can2020_trt <- lm(climateworry ~ climexp_trt, data=can2020)
coeftest(ate_worry_can2020_trt, vcov=vcovHC(ate_worry_can2020_trt, type="HC1"))


### USA 2020 ##########

# Read in US 2020 survey data
usa2020 <- read.csv("usdata_summer2020.csv")

# Code treatment groups 
usa2020$climexp_cont <- ifelse(usa2020$climateexp=="A", 1, 0)   
usa2020$climexp_t_txt <- ifelse(usa2020$climateexp=="B", 1, 
                                ifelse(usa2020$climateexp=="A", 0, NA))
usa2020$climexp_t_img <- ifelse(usa2020$climateexp=="C", 1, 
                                ifelse(usa2020$climateexp=="A", 0, NA))
usa2020$climexp_trt <- ifelse(usa2020$climateexp=="B" | usa2020$climateexp=="C", 1, 0)

# DV - in the US data we have both climate worry and climate priority in 2020
usa2020$climate_worry <- recode(usa2020$climate_worry, "Not at all worried"=0, "Not very worried"=1, "Somewhat worried"=2, "Very worried"=3)
usa2020$climate_priority <- recode(usa2020$climate_priority, "Low"=0, "Medium"=1, "High"=2, "Very high"=3)

#ATEs on 2020 US climate worry 
ate_worry_usa2020_txt <- lm(climate_worry ~ climexp_t_txt, data=usa2020)
coeftest(ate_worry_usa2020_txt, vcov=vcovHC(ate_worry_usa2020_txt, type="HC1"))
ate_worry_usa2020_img <- lm(climate_worry ~ climexp_t_img, data=usa2020)
coeftest(ate_worry_usa2020_img, vcov=vcovHC(ate_worry_usa2020_img, type="HC1"))
ate_worry_usa2020_trt <- lm(climate_worry ~ climexp_trt, data=usa2020)
coeftest(ate_worry_usa2020_trt, vcov=vcovHC(ate_worry_usa2020_trt, type="HC1"))

#ATEs on 2020 US climate priority
ate_priority_usa2020_txt <- lm(climate_priority ~ climexp_t_txt, data=usa2020)
coeftest(ate_priority_usa2020_txt, vcov=vcovHC(ate_priority_usa2020_txt, type="HC1"))
ate_priority_usa2020_img <- lm(climate_priority ~ climexp_t_img, data=usa2020)
coeftest(ate_priority_usa2020_img, vcov=vcovHC(ate_priority_usa2020_img, type="HC1"))
ate_priority_usa2020_trt <- lm(climate_priority ~ climexp_trt, data=usa2020)
coeftest(ate_priority_usa2020_trt, vcov=vcovHC(ate_priority_usa2020_trt, type="HC1"))


### FIGURE 4 IN PAPER ##########

#Extract coeftest results
plot_canworry_2020 <- broom::tidy(coeftest(ate_worry_can2020_trt, vcov=vcovHC(ate_worry_can2020_trt, type="HC1")), conf.int=T)
plot_usaworry_2020 <- broom::tidy(coeftest(ate_worry_usa2020_trt, vcov=vcovHC(ate_worry_usa2020_trt, type="HC1")), conf.int=T)
plot_canworry_2021 <- broom::tidy(coeftest(ate_worry_can_trt, vcov=vcovHC(ate_worry_can_trt, type="HC1")), conf.int=T)
plot_usaworry_2021 <- broom::tidy(coeftest(ate_worry_usa_trt, vcov=vcovHC(ate_worry_usa_trt, type="HC1")), conf.int=T)
plot_canpriority_2021 <- broom::tidy(coeftest(ate_priority_can_trt, vcov=vcovHC(ate_priority_can_trt, type="HC1")), conf.int=T)
plot_usapriority_2021 <- broom::tidy(coeftest(ate_priority_usa_trt, vcov=vcovHC(ate_priority_usa_trt, type="HC1")), conf.int=T)

#Add DV labels with sample
plot_canworry_2020[2,1] <- "Climate Worry (2020)"
plot_canworry_2020 <- plot_canworry_2020[2, c(1:2, 6:7)]
plot_usaworry_2020[2,1] <- "Climate Worry (2020)"
plot_usaworry_2020 <- plot_usaworry_2020[2, c(1:2, 6:7)]
plot_canworry_2021[2,1] <- "Climate Worry (2021)"
plot_canworry_2021 <- plot_canworry_2021[2, c(1:2, 6:7)]
plot_usaworry_2021[2,1] <- "Climate Worry (2021)"
plot_usaworry_2021 <- plot_usaworry_2021[2, c(1:2, 6:7)]
plot_canpriority_2021[2,1] <- "Climate Priority (2021)"
plot_canpriority_2021 <- plot_canpriority_2021[2, c(1:2, 6:7)]
plot_usapriority_2021[2,1] <- "Climate Priority (2021)"
plot_usapriority_2021 <- plot_usapriority_2021[2, c(1:2, 6:7)]

#Make these into a dataframe for plotting
plotdata_svyexp_coefs <- rbind(plot_canworry_2020, plot_usaworry_2020, plot_canworry_2021, plot_usaworry_2021, plot_canpriority_2021, plot_usapriority_2021)

#Add country column for faceting & adjust factor order
plotdata_svyexp_coefs$country <- c("Canada", "United States", "Canada", "United States", "Canada", "United States")
plotdata_svyexp_coefs$country <- factor(plotdata_svyexp_coefs$country, levels=c("United States", "Canada"))

#Create plot
svyexp_coefplot_facet <- ggplot(data=plotdata_svyexp_coefs, 
                                mapping=aes(x=term, y=estimate)) +
  geom_point(size=1.5) +
  ylab("Average treatment effect") +
  xlab("Dependent variable") +
  facet_grid(cols=vars(country)) +
  ggtitle("Exponential growth vignette experiment results") +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.1) +
  geom_hline(yintercept=0, linetype="dashed", color="darkgray") +
  coord_flip() +
  theme_bw()
#Save plot
ggsave(svyexp_coefplot_facet, file=paste0(figfolder,"svyexp_coefs_all_facet.png"), width=8, height=4, units="in",dpi=600)


### TABLE 6 IN THE SI ##########

# Get robust standard errors:
worry_can2020_cov <- vcovHC(ate_worry_can2020_trt, type="HC1")
worry_can2020_se <- sqrt(diag(worry_can2020_cov))
worry_can2021_cov <- vcovHC(ate_worry_can_trt, type="HC1")
worry_can2021_se <- sqrt(diag(worry_can2021_cov))
worry_usa2020_cov <- vcovHC(ate_worry_usa2020_trt, type="HC1")
worry_usa2020_se <- sqrt(diag(worry_usa2020_cov))
worry_usa2021_cov <- vcovHC(ate_worry_usa_trt, type="HC1")
worry_usa2021_se <- sqrt(diag(worry_usa2021_cov))

stargazer(ate_worry_can2020_trt, ate_worry_can_trt, ate_worry_usa2020_trt, ate_worry_usa_trt,
          title = "Average treatment effects of analogizing climate change to the exponential growth of COVID-19 on respondents’ level of worry about climate change.",
          label = "tab:svyexp_worry_bundled",
          covariate.labels = c("Average treatment effect"),
          model.names = FALSE,
          multicolumn = FALSE,
          dep.var.caption  = "DV: worry about climate change",
          dep.var.labels.include = FALSE,
          column.labels = c("CAN: May 2020", "CAN: Jun. 2021", "US: May 2020", "US: Feb. 2021"),
          se = list(worry_can2020_se, worry_can2021_se, worry_usa2020_se, worry_usa2021_se),
          add.lines = list(c("Observations", "892", "1,040", "1,049", "1,609")),
          omit.stat = c("all"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = "Robust standard errors in parentheses.",
          digits = 3,
          header = FALSE,
          out = paste0(tabfolder,"svyexp_worry_bundled.tex"))

### TABLE 7 IN THE SI ##########

# Get robust standard errors:
priority_can2021_cov <- vcovHC(ate_priority_can_trt, type="HC1")
priority_can2021_se <- sqrt(diag(priority_can2021_cov))
priority_usa2021_cov <- vcovHC(ate_priority_usa_trt, type="HC1")
priority_usa2021_se <- sqrt(diag(priority_usa2021_cov))

stargazer(ate_priority_can_trt, ate_priority_usa_trt,
          title = "Average treatment effects of analogizing climate change to the exponential growth of COVID-19 on respondents’ belief that government should prioritize climate change.",
          label = "tab:svyexp_priority_bundled",
          covariate.labels = c("Average treatment effect"),
          model.names = FALSE,
          multicolumn = FALSE,
          dep.var.caption  = "DV: prioritization of climate change",
          dep.var.labels.include = FALSE,
          column.labels = c("CAN: Jun. 2021", "US: Feb. 2021"),
          se = list(priority_can2021_se, priority_usa2021_se),
          omit.stat = c("all"),
          add.lines = list(c("Observations", "1,040", "1,611")),
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = "Robust standard errors in parentheses.",
          digits = 3,
          header = FALSE,
          out = paste0(tabfolder,"svyexp_priority_bundled.tex"))

### TABLE 8 IN THE SI ##########

#### Generate full model with unbundled treatments

#USA 2021
usasvy$climexptxt <- ifelse(usasvy$climateexperiment=="B", 1, 0) #Text only treatment
usasvy$climexpimg <- ifelse(usasvy$climateexperiment=="C", 1, 0) #image treatment
ate_worry_usa_full <- lm_robust(climateworry~climexptxt+climexpimg, data=usasvy, se_type="HC1")
summary(ate_worry_usa_full)
#Test whether the treatments are different
linearHypothesis(ate_worry_usa_full,hypothesis.matrix=c("climexptxt=climexpimg"),test="F")

#Canada 2021
cansvy$climexptxt <- ifelse(cansvy$climateexperiment=="B", 1, 0) #Txt only trtmt
cansvy$climexpimg <- ifelse(cansvy$climateexperiment=="C", 1, 0) #Img trtmt
ate_worry_can_full <- lm_robust(climateworry~climexptxt+climexpimg, data=cansvy, se_type="HC1")
summary(ate_worry_can_full)
#Test whether the treatments are different
linearHypothesis(ate_worry_can_full,hypothesis.matrix=c("climexptxt=climexpimg"),test="F")

#USA 2020
usa2020$climexptxt <- ifelse(usa2020$climateexp=="B", 1, 0) #Txt only trtmt
usa2020$climexpimg <- ifelse(usa2020$climateexp=="C", 1, 0) #Img trtmt
ate_worry_usa2020_full <- lm_robust(climate_worry~climexptxt+climexpimg, data=usa2020, se_type="HC1")
summary(ate_worry_usa2020_full)
#Test whether treatments are different
linearHypothesis(ate_worry_usa2020_full,hypothesis.matrix=c("climexptxt=climexpimg"),test="F")

#Canada 2020
can2020$climexptxt <- ifelse(can2020$GROUP_EXP=="Group 2", 1, 0) #Txt only trtmt
can2020$climexpimg <- ifelse(can2020$GROUP_EXP=="Group 3", 1, 0) #Img trtmt
ate_worry_can2020_full <- lm_robust(climateworry~climexptxt+climexpimg, data=can2020, se_type="HC1")
summary(ate_worry_can2020_full)
#Test whether treatments are different
linearHypothesis(ate_worry_can2020_full,hypothesis.matrix=c("climexptxt=climexpimg"),test="F")

#Only survey waves with significant effects: CAN2020 (when bundled) and US2021
#Generate table:
texreg(list(ate_worry_can2020_full, ate_worry_usa_full),
       include.ci=FALSE,
       custom.model.names=c("CAN: May 2020","US: Feb. 2021"),
       custom.coef.names=c("Control Group Mean", "Text only trt. effect","Text + image trt. effect"),
       caption="Average treatment effects of analogizing climate change to the exponential growth of COVID-19 on respondents' level of worry about climate change. ATEs are presented for both the text-only treatment and a treatment that included images added to the text showing graphical representations of exponential growth. We present results in this table for the survey waves in which we find significant treatment effects.",
       file=paste0(tabfolder,"svyexp_worry_unbundled.tex"),
       label="tab:svyexpworry_unbundled",
       caption.above=TRUE,
       float.pos="h")

### TABLE 9 IN SI ##########

#USA 2021
ate_priority_usa_full <- lm_robust(climatepriority~climexptxt+climexpimg, data=usasvy, se_type="HC1")
summary(ate_priority_usa_full)
#Test whether treatments are different
linearHypothesis(ate_priority_usa_full,hypothesis.matrix=c("climexptxt=climexpimg"),test="F")

#Canada 2021
ate_priority_can_full <- lm_robust(climatepriority~climexptxt+climexpimg, data=cansvy, se_type="HC1")
summary(ate_priority_can_full)
#Test whether treatments are different
linearHypothesis(ate_priority_can_full,hypothesis.matrix=c("climexptxt=climexpimg"),test="F")

#USA 2020
ate_priority_usa2020_full <- lm_robust(climate_priority~climexptxt+climexpimg, data=usa2020, se_type="HC1")
summary(ate_priority_usa2020_full)
#Test whether treatments are different
linearHypothesis(ate_priority_usa2020_full,hypothesis.matrix=c("climexptxt=climexpimg"),test="F")

#Generate table:
texreg(list(ate_priority_usa_full),
       include.ci=FALSE,
       custom.model.names=c("US Jun. 2021"),
       custom.coef.names=c("Control group mean", "Text only trt. effect","Text + image trt. effect"),
       caption="Average treatment effects of analogizing climate change to the exponential growth of COVID-19 on respondents' belief that government should make climate change a priority. ATEs are presented for both the text-only treatment and a treatment that included images added to the text showing graphical representations of exponential growth. We present results in this table for the survey wave in which we found significant effects.",
       file=paste0(tabfolder,"svyexp_priority_unbundled.tex"),
       label="tab:svyexppriority_unbundled",
       caption.above=TRUE,
       float.pos="h")
### Mini-meta analysis #### 
#Construct dataframe to feed into `metafor`

study <- c("Canada (May 2020)", "Canada (Jun. 2021)", "US (May 2020)", "US (Feb. 2021")
#N in treatment group
n_treat <- c(sum(can2020$climexp_trt[!is.na(can2020$climateworry)], na.rm=T),
             sum(cansvy$climexp_trt[!is.na(cansvy$climateworry)], na.rm=T),
             sum(usa2020$climexp_trt[!is.na(usa2020$climate_worry)], na.rm=T),
             sum(usasvy$climexp_trt[!is.na(usasvy$climateworry)], na.rm=T))
#N in control group
n_cont <- c(sum(can2020$climexp_cont[!is.na(can2020$climateworry)], na.rm=T),
            sum(cansvy$climexp_cont[!is.na(cansvy$climateworry)], na.rm=T),
            sum(usa2020$climexp_cont[!is.na(usa2020$climate_worry)], na.rm=T),
            sum(usasvy$climexp_cont[!is.na(usasvy$climateworry)], na.rm=T))
#Mean worry in treated group
mean_treat <- c(mean(can2020$climateworry[can2020$climexp_trt==1], na.rm=T),
                mean(cansvy$climateworry[cansvy$climexp_trt==1], na.rm=T),
                mean(usa2020$climate_worry[usa2020$climexp_trt==1], na.rm=T),
                mean(usasvy$climateworry[usasvy$climexp_trt==1], na.rm=T))
#Mean worry in control group
mean_cont <- c(mean(can2020$climateworry[can2020$climexp_cont==1], na.rm=T),
               mean(cansvy$climateworry[cansvy$climexp_cont==1], na.rm=T),
               mean(usa2020$climate_worry[usa2020$climexp_cont==1], na.rm=T),
               mean(usasvy$climateworry[usasvy$climexp_cont==1], na.rm=T))
#SD treatment group
sd_treat <- c(sd(can2020$climateworry[can2020$climexp_trt==1], na.rm=T),
              sd(cansvy$climateworry[cansvy$climexp_trt==1], na.rm=T),
              sd(usa2020$climate_worry[usa2020$climexp_trt==1], na.rm=T),
              sd(usasvy$climateworry[usasvy$climexp_trt==1], na.rm=T))
#SD control group
sd_cont <- c(sd(can2020$climateworry[can2020$climexp_cont==1], na.rm=T),
             sd(cansvy$climateworry[cansvy$climexp_cont==1], na.rm=T),
             sd(usa2020$climate_worry[usa2020$climexp_cont==1], na.rm=T),
             sd(usasvy$climateworry[usasvy$climexp_cont==1], na.rm=T))
#Combine into dataframe
metadata <- data.frame(study, n_cont, n_treat, mean_cont, mean_treat, sd_cont, sd_treat)

#Apply metafor scaling function for standardized mean difference 
metadata <- escalc(measure="SMD", data=metadata,
                   n1i = n_treat, n2i = n_cont,
                   m1i = mean_treat, m2i = mean_cont,
                   sd1i = sd_treat, sd2i = sd_cont,
                   append = TRUE)

#Calculate random effects model to get overall effect
minimeta_mod <- rma(yi, vi, data=metadata)
summary(minimeta_mod)

#Generate Forest plot
png(paste0(figfolder,"minimeta_forestplot.png"),
    width=8, height=4, units="in", res=600)
forest(minimeta_mod, slab=metadata$study, header=c("Study", "Std. Mean Diff. [95% CI]"), digits=3L)
dev.off()
